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We use a Id model of a superfluid based on the Gross-Pitaevskii Lagrangian to illustrate a general 
numerical method designed to find quantum tunneling rates in extended bosonic systems. Specif- 
ically, we study fiow past an obstacle and directly solve the imaginary time dynamics to find the 
"bounce" solution connected with the decay of the metastable laminar state via phase slip nucle- 
ation. The action for the tunneling configuration goes to zero at the threshold (in superfluid velocity) 
for classical production of these slips. Applications to other processes are briefly discussed. 

PACS numbers: 64.60. Qb, 67.40.-w, 67.40. Vs 
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The subject of quantum tunneling in extended systems 
continues to be of great interest in a variety of condensed 
matter applications. Examples include superflows (either 
superfluids or superconductors) in which dissipation can 
be induced by the creation of phase slips via tunneling. 
Often, the rates of these processes are calculated by fo- 
cusing on one relevant degree of freedom (perhaps the 
center-point of a specific ansatz for the quantum field) 
and approximately integrating out the remainder of the 
degrees of freedom; this approach thereby reduces the 
problem to quantum mechanics of a single particle. This 
reduction precludes one from understanding the extent 
to which re-arrangement of the (usually nonlinear) field 
degrees of freedom can lead to more advantageous tun- 
neling paths. 

In this work, we offer a more direct approach to 
the evaluation of quantum tunneling rates in extended 
bosonic systems. Specifically, we use the idea popular- 
ized by Coleman and Callan that quantum tunneling can 
be studied by finding a specific solution to the classical 
equations of motion continued to imaginary time, the so- 
called "bounce". To illustrate our method, we explicitly 
study the nucleation of phase slips in a Id model of a 
superfluid flowing past an obstacle. We show how the 
bounce solution can be constructed, how it determines 
the nucleated state of the entire field and how one can 
compute the action for this solution, thereby finding the 
semiclassical rate. At the end, we discuss extensions to 
other bosonic systems as well as comment on the modi- 
fications that will be required to study fermions. 

To model a superflow with velocity v past a Id obstacle 
we use the Gross-Pitaevskii Lagrangian density j|] 



C = ihij)dtip-— \dxip\^- 
zm 



The first three terms arise directly from considering 
bosons interacting with each other via S function poten- 
tials of strength g, in the presence of an obstacle treated 
as a position-dependent chemical potential fJ-{x). The 



last piece comes from working in a frame in which the 
obstacle is at rest. We will explicitly choose /x(x) to 
have the form of a square well of width A and depth 
rjfiQ (rj < 1), fiQ being the uniform chemical potential 
outside; this will give rise to a depletion of the conden- 
sate density in the well. In this model, the natural unit 
of length is the coherence length ^ = h/y/rnjlo, of veloc- 
ity is the sound velocity, c = ^ IJ-o/m, and of condensate 
amplitude is ^/po = y/io/y. A dimensionless param- 
eter, a = poS,'^ {d = 1 here), analogous to 1/h, then 
determines the importance of quantum fluctuations; we 
assume a ^ 1. The classical dynamics derived from 
the above Lagrangian is the Galilean-transformed non- 
linear Schrodinger equation (NLSE). In natural units, 
with {((x) = p{x)/fio), 



In this equation, ip is often interpreted as represent- 
ing "the wavefunction of the condensate". The spatial 
boundary condition is taken to be one where the ip tends 
to a unimodular constant far away from the obstacle. 

In recent work by Frisch et al |^], the NLSE in two 
spatial dimensions was utilized to study the problem of 
2d superflow past a circular obstacle. They found that 
there exists a critical velocity such that for v < Vc, one 
is able to find time-independent solutions correspond- 
ing to a laminar flow around the obstacle, whereas for 
V > Vc only time-dependent solutions exist. The lat- 
ter solutions correspond to a superflow that periodically 
nucleates vortex-antivortex pairs in the proximity of the 
disk; these pairs are subsequently advected by the flow. 
Similar results were found in Ref. [| for a superflow in 
a 2d channel containing a constriction. The Id analog 
of these vortex-antivortex pairs are phase slips. These 
are events localized in both space and time where the 
amplitude of the condensate goes to zero and the phase 
slips by 27r. Accordingly, we investigated numerically the 
existence of time-independent solutions of this equation; 
this was done by writing the time-independent NLSE in 
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a Id grid and using Newton's method complemented by 
Powell's hybrid method |^,|| to solve the resulting sys- 
tem of nonlinear algebraic equations. In direct analogy 
to the 2d findings, we found that, associated with each 
well depth and width, there is a critical value of the su- 
perflow velocity, Vc, above which such solutions do not 
exist. 

The above result implies that the solution of the NLSE 
above is time-dependent. To find the nature of this 
solution, the time-dependent NLSE was integrated using 
an operator splitting method |^ accurate to 0{dt) and 
unconditionally stable. An example of what we found 
is shown in Fig. 1. We observed the nucleation of the 
phase slip, typically at the edge of the well where \dxC\ is 
greatest, and the subsequent emission of a soliton. This 
soliton, which propagates outside the well, is a member 
of a family of solutions of the NLSE with C = 1 0- Up 
to an overall phase, 



ipui^, t) = VT— u^tanh \J 1 — [a; — (u + v)t 



-|- iu. 



This is a "dark" soliton: the amplitude dips to \u\ at its 
center. It carries a momentum P = —2u\/l — and 
accounts for a phase slip of tt — 2sin^^ u (see Fig. 1). 




FIG. 1. The time evolution for v > Vc, showing amplitude 
(solid line) and phase (dashed line) of i/" at 7 times between 
and 24 ^/c. The phase slip occurs at t « 12 ^/c, pre- 
cisely at the right edge of the square well. Afterwards, a 
soliton is propagates to the right. The well width was 1.6^ 
and the chemical potential inside was 0.3/^o. The velocity was 
V = 0.4 c > Wc = 0.3 c. was fixed at the ends. Inset: critical 
velocity versus well depth for wells of width 3.2^ (squares) 
and 1.6^ (triangles). 

The existence of time-independent solutions of the 
NLSE only for v < Vc suggests a picture where the 



parameter v creates a barrier that prevents the classi- 
cal decay of the metastable, time-independent state 
This barrier disappears when v > Vc and thus only time- 
dependent solutions can be found. However, the steady 
flow state is only metastable and therefore can decay via 
quantum fluctuations. This is similar to the "fate of the 
false vacuum" problem studied by Coleman and Callan 
One can therefore use their formalism which relates 
the the tunneling rate (to leading order in a semiclassical 
expansion) to the classical action of the periodic "bounce 
solution" , which proceeds in imaginary time from the 
metastable state to the nucleated field and back again. 
With the exception of one paper which studies nuclear 
fission jl^], the possibility that one could directly solve 
for these bounce solutions does not seem to have been 
previously realized. 

In detail, we consider a coherent state path integral 
representation of the matrix element of the evolution op- 
erator for imaginary time. 



(V'ole 



--HT/h 



1^0 



V[^,i^] exp[-^E[^,Vi]/?i] 



V'(-r/2)=^o 



Here, |-0o) is the metastable coherent state wavefunction. 
The Euclidean action is {L is the spatial size of the sys- 
tem) 

T/2 L/2 

S^[i^, ^] = Po ^ Jdrjdx [i;dr^ + E{^, V)] , (1) 

-T/2 -L/2 

E{^, ^P) = i(a,^)(9,i^) + - C(^)]' - ^v^d,^ . (2) 

The bounce solution is a classical configuration which 
extremizes this action; we easily obtain the bounce equa- 
tions for the fields ip and 'ip: 



dri> = - ^dl^p + [Tpip - C,[x)]ip - ivdxtp 
drip = — hdl'ip + [tpip — Ci^)]''P + ii'dxip 



(3) 
(4) 



with to ip{-T/2) = -(/lo and i'{T/2) = li:^. The spatial 
boundary conditions are the same as before. These equa- 
tions imply V'(t) = ■0(— r)*, hence iP{t) = 'iP*{t) - the 
density %pip is real - only at three time slices: r = i^T, 
when ip = ipQ, and r = 0, when the field emerges from the 
barrier and the phase slip has been nucleated |11|. It is 
crucial to recognize that we do not need to make any as- 
sumption regarding the nature of the nucleated state, as 
the solutions of the equations will directly and uniquely 
determine this field configuration. 

These equations can not be solved by integrating for- 
ward in time. This is because the ip field has negative 
diffusivity and hence forward in time integration is ill- 
posed; alternatively, eliminating ip using the above rela- 
tion will give an explicitly non-local (in time) system. 
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FIG. 2. Contour plot of the amplitude (solid) and phase 
(dashed) of ip{x,T) for the bounce solution. Here v = 0.1c 
while Vc — 0.3 c. Heavy lines arc where the phase crosses —2tt 
and do not correspond to any singularity. The well parameters 
are the same as in Fig. 1. 

Even were this not to be the case, the fact that the 
bounce is a saddle-point as opposed to a minimum ren- 
ders time-integration difficult. Instead, we solved those 
equations iteratively by discretizing them on a space-time 
lattice and using the Newton-Powell method mentioned 
above to find the zero of the resulting system of nonlinear 
algebraic equations. This iterative scheme does not care 
about time-locality, nor does it assume stability. 

In solving the bounce equations on a space-time lat- 
tice one must choose a time interval T long enough to 
include all the structure around r = 0, as it is that part 
of the bounce that contributes the most to the action 
and that gives the nucleated state. We set up the equa- 
tions in such a way that at r = —T/2 we had the exact 
metastable state, while at once we also enforced the re- 
lation V'(t) = %//*{— t). Because T is finite, the solution 
obtained is not exactly periodic - tp at t = T/2 is slightly 
different from tp at t = —T/2 - but the difference can 
be made arbitrarily small by increasing T. The energy 
of the solution, however, was nearly constant and equal 
to the metastable state energy, as is the case for the true 
bounce solution. We also work with a finite spatial size. 
Whereas in the infinite system size limit we can have both 
ip — > e*^" and dxip = as a; ^ ±oo, in our case we had 
to choose between these two options. We found that for 
smaller system sizes the gradient condition worked best. 
Typically we had T ^ 34 ^/c and L ^ 20^ in a grid of 
100 spatial points by 50 time points. This resulted in 
a systcim of 10,000 real nonlinear algebraic equations for 
the real and imaginary parts of ijj. The iterative scheme 
converges rapidly to a well-behaved solution when an ini- 



FIG. 3. The surface plot of the amplitude of the field ip of 
Fig. 2. The depletion of the condensate inside the well and 
in the vortex cores is apparent. 

tial guess with the correct topology is used to start the 
process. 

As expected, we found a solution that precisely corre- 
sponds to the quantum nucleation of a phase slip from 
the time-independent state of the given well at superflow 
velocity {v < Vc)- The field ip{x,T) corresponding to a 
phase slipped state at r = has a singular point in nega- 
tive time around which its phase winds counter-clockwise 
and another at positive time with clockwise phase wind- 
ing. This vortex-antivortex pair in space-time results in 
a state at t = with a phase slip of 2t: with respect to 
the metastable state at r = ±T/2. The amplitude and 
phase of a typical solution tp are shown in Figs. 2 and 3. 
Note that the phase singularities are not symmetrically 
displaced about t = - ip and ip are only related by 
complex conjugation for t — 0, ±T/2. Precisely at r = 
the configuration is interpretable as one for the ordinary 
{i.e. real-time) NLSE and hence can be treated as the 
initial condition for the subsequent real-time evolution of 
the nucleated phase slip. When this is done, the result 
is quite similar to what we have already seen occurs for 
the classically nucleated slip at u > Vc', that is, a soli- 
ton is formed which propagates to the right, transferring 
energy from the obstacle region towards infinity. 

The decay rate of the metastable state can be found 

using r c>c exp[— AS'e/Ti], where A^e is the difference be- 
tween the bounce action and the metastable state action. 
Since the dynamics in imaginary time conserves energy, 
we have 



r (X exp — po ^ dr dx tpdrip 



(5) 
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We solved the bounce equations for various superflow ve- 
locities with a given well and computed T according to 
equation ^. The behavior is depicted in Fig. 4; note 
that the action approaches zero as we reach the classical 
threshold. As v approaches Vc the vortex-antivortex pair 
separation decreases until eventually no time-dependent 
solution of equations || can be found for v > Vc- 



-2 - 



o 



-6 - 



-8 - 



-10 



-12 



0.2 



0.4 



0.6 



0.8 



FIG. 4. An example of the dependence of the semiclassical 
decay rate on the superflow velocity. The well is the same as 
Figs. 2 and 3. For v > Vc no metastable state exists. 

Although we have described our method for a specific 
example, it should be clear that the same approach can 
be used to compute tunneling rates and nucleated states 
for any bosonic system. Applications we have in mind 
include higher dimensional superfluid-obstacle problems 
(such as the quantum analog of the 2d classical nucle- 
ation of vortex-antivortex pairs discussed above) and the 
decay of the metastable Bose-Einstein condensate state 
for trapped alkali atoms with attractive interactions [ p2| . 
One can also imagine straightforward extensions to the 
traditional field of macroscopic quantum tunneling where 
one could study the tunneling of, say, a two level sys- 
tem coupled to a nonlinear set of continuous degrees of 
freedom. Finally, extensions to fermionic systems (either 
charged or uncharged) are more challenging. Integrating 
out the fermions typically leaves one with a bosonic ac- 
tion in the form of some functional determinant and the 
resultant equation of motion involves successfully eval- 
uating the inverse of a complicated operator. Whether 
it is feasible to directly find the bounce solution for this 
case (with no further approximations) is currently under 
investigation. 
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